function [results, standard_deviation, img_points] = read_rect_target(filepath, gridsize, varargin)
% READ_RECT_TARGET, a MATLAB target measuring script.
% Copyright (C) 2012 Glenn D. Sweeney
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.
% 
%
%
%
% Description:
% READ_RECT_TARGET is a function designed for measuring average digital
% counts for each patch in an image of a rectangular target. It uses user
% input to locate the target in the image, and is capable of processing
% oblique images of targets.
%
% Usage:
% After calling READ_RECT_TARGET, an image of the test target will appear.
% The user then locates the target by clicking on the four corners of the
% rectangular grid, in any order.
% 
% Syntax:
% READ_RECT_TARGET(filepath, [num_rows num_cols])
% READ_RECT_TARGET(..., param1, val1, param2, val2,...)
% output = READ_RECT_TARGET(...)
%
% Mandatory Arguments:
% 
%   filepath            This argument specifies the location of the image
%                       (or folder of images, see below) to be processed
%
%   num_rows            The number of rows in the rectangular test target
%
%   num_cols            The number of columns in the rectangular test
%                       target
%
% Optional Arguments:
%   Optional arguments are supplied using a 'parameter', (value) syntax,
%   where function parameters may be specified in any order, and then the
%   paired value (if any, depending on the parameter) is passed in as the
%   next function argument. For example, if the parameter 'CropFactor' is
%   to be set to 0.5, the calling syntax would be:
%
%   READ_RECT_TARGET(filepath, [num_rows, num_cols], 'CropFactor', 0.5)
%
%   A table of parametrs with a description of acceptable values is
%   provided below:
%
% _________________________________________________________________________
% |                    |                                                  |
% |    'CropFactor'    | The value immediately following 'CropFactor'     |
% |                    | should be a decimal between zero and one. This   |
% |                    | value denotes the proportion  of each patch to   |
% |                    | be removed from averaging.                       |
% |                    | Default value is 0.25.                           |
% |                    | Usage: read_rect_target(..., 'CropFactor', 0.25) |
% |                    |                                                  |
% |     'Mosaic'       | Mosaic allows a user to input a 2x2 integer      |
% |                    | array denoting the mosaic pattern present in a   |
% |                    | grayscale image. The pattern must consist only   |
% |                    | of 1, 2, and 3, corresponding to R,G,B channels. |
% |                    | An example is as follows:                        |
% |                    |      bayerarr =    ['2' '1'; ...                 |
% |                    |                     '3' '2']                     |
% |                    |      read_rect_target(...,'Mosaic', bayerarr)    |
% |                    | If 'Mosaic' is not specified, it is assumed that |
% |                    | the image has already been demosaiced.           |
% |                    |                                                  |
% |    'Normalize'     | Normalize scales function output from 0 to 1 for |
% |                    | display purposes. It does not affect output data.|
% |                    | This tag is NOT followed by a value.             |
% |                    |                                                  |
% |     'Folder'       | Specifying the folder parameter will process a   |
% |                    | folder of identically-composed images. This tag  |
% |                    | is not followed by a value.                      |
% |                    | Usage: read_rect_target(folderpath, [r,c], ...,  |
% |                    |          'Folder')                               |
% |                    |                                                  |
% |    'ImgPoints'     | Allows the user to pass in a set of predetermined|
% |                    | control points for processing (same format as    |
% |                    | imgpoints return value)                          |
% |                    |                                                  |
% |  'SuppressDraw'    | When this flag is set, no figures will be created|
% |                    | for output .This parameter is not followed by a  |
% |                    | value.                                           |
% |____________________|__________________________________________________|
%
% Output Values:
%   READ_RECT_TARGET returns three sets of values. The first is an RxC
%   array, where each array location has a vector of average digital counts
%   for each channel. Accordingly, for an RGB image, this output will be a
%   3d array, indexable by: (row, col, chan).
%   Results will always be doubles, but scaled to match the bit depth of
%   the input image.
%   Usage: result = READ_RECT_TARGET(...)
%
%   Optionally, READ_RECT_TARGET will return standard deviations for each
%   patch, in the same format as patch digital counts.
%   Usage: [digital_count, stdev, ~] = READ_RECT_TARGET(...)
%
%   Finally, if many images are going to be processed with the same target,
%   but for some reason the 'Folder' parameter is unsatisfactory, the
%   selected coordinates are also returned, and can be passed into a
%   subsequent call of READ_RECT_TARGET through the 'ImgPoints' parameter.
%   Usage: [digital_count, stdev, impoints] = READ_RECT_TARGET(...)

% Default variable states
crop_factor = .25;
mosaic_array = 0;
peak_normalize = 0;
folder_flag = 0;
points_flag = 0;
draw_flag = 1;


% Input argument parsing
i = 1;
while i <= numel(varargin)
    
    %Check the argument to see if it matches a parameter. If so, set the
    %next values accordingly
    switch varargin{i}
        case 'CropFactor'
            crop_factor = varargin{i+1};
        case 'Mosaic'
            mosaic_array = varargin{i+1};
        case 'Normalize'
            peak_normalize = 1;
            i = i-1;
        case 'Folder'
            folder_flag = 1;
            i = i-1;
        case 'ImgPoints'
            points_flag = 1;
            img_points = varargin{i+1};
        case 'SuppressDraw'
            draw_flag = 0;
            i = i-1;
        otherwise
            error('Unknown function argument.');
    end
    i = i+2;
end

% Get the number of rows and columns
num_rows = gridsize(1);
num_cols = gridsize(2);

% Go through the list of files (if filepath is a single image, this will
% only be populated by a single file)
file_list = dir(filepath);

% On a mac, hidden files begin with '.'. Ignore these files.
while strcmp(file_list(1).name(1), '.')
    file_list = file_list(2:end);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Find Target Corners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Read in the target image (for point selection ONLY)
if folder_flag
    input = imread([filepath file_list(1).name]);
else
    input = imread(filepath);
end

% Process the mosaic by splitting across color channels. Populate all other
% pixels with -1 to avoid averaging them in.
if mosaic_array == 0
    target = input;
else
    tar_size = size(input);
    tar_height = tar_size(1);
    tar_width = tar_size(2);
    target = zeros([tar_size, 3])-1;
    target(1:2:tar_height, 1:2:tar_width, mosaic_array(1,1)) = input(1:2:tar_height, 1:2:tar_width);
    target(2:2:tar_height, 2:2:tar_width, mosaic_array(2,2)) = input(2:2:tar_height, 2:2:tar_width);
    target(1:2:tar_height, 2:2:tar_width, mosaic_array(1,2)) = input(1:2:tar_height, 2:2:tar_width);
    target(2:2:tar_height, 1:2:tar_width, mosaic_array(2,1)) = input(2:2:tar_height, 1:2:tar_width);
end



% Calculate the size of the target
tar_size = size(target);
tar_height = tar_size(1);
tar_width = tar_size(2);

% Test for number of color channels
if numel(tar_size) == 3
    num_channels = tar_size(3);
else
    num_channels = 1;
end

% Define the function output
output = zeros(num_rows, num_cols, num_channels);
stdev = zeros(num_rows, num_cols, num_channels);

results = zeros([num_rows, num_cols, num_channels, numel(file_list)]);
standard_deviation = zeros([num_rows, num_cols, num_channels, numel(file_list)]);

% If corners were not specified by the user, find them now.
if ~points_flag

    % Display the image, and read in 4 corners specified by user
    if peak_normalize == 0
        imshow(target);
    else
        imshow(target./max(target(:)));
    end
    hold on
    num_good_points = 0;

    % first row: x   Second row: y    3rd row: homogeneous coordinate
    img_points = zeros(3, 4);

    % Size of the corner dots to be used (ceil to ensure at least one pixel)
    % Defined as 1/100 the smaller dimension
    dot_size = ceil(min([tar_height tar_width])/100);

    % Get target registration points
    while num_good_points < 4

        [x_click y_click button] = ginput(1);
        % If it was left mouse click:
        if button == 1
            % Show the pretty red squares
            fill([x_click+dot_size x_click+dot_size x_click-dot_size x_click-dot_size], ...
                [y_click+dot_size y_click-dot_size y_click-dot_size y_click+dot_size], 'r');
            num_good_points = num_good_points + 1;
            img_points(1, num_good_points) = x_click;
            img_points(2, num_good_points) = y_click;
            img_points(3, num_good_points) = 1; % Generating homogeneous coordinates
            drawnow;
        end
    end

    close all;

end

% Sorts the four image points so they go clockwise from top left
img_points = sortrows(img_points', 2)';
img_points(:, 1:2) = sortrows(img_points(:, 1:2)', 1)';
img_points(:, 3:4) = sortrows(img_points(:, 3:4)', -1)';

% Target points assume upper left corner is zero, and bottom right is
% (num_rows, num_cols), sorted clockwise from top left
target_points = [0 num_cols num_cols    0; ...
                 0    0     num_rows num_rows; ...
                 1    1        1        1];

% Generate A matrix for each correspondance:
% Projecting from the src (source) system into the dst (destination) system
A = [];
for i = 1:4
    src = target_points(:,i);
    dst = img_points(:,i);
    
    A = vertcat( A, [0,0,0,   -1*src',   dst(2)*(src');...
        src',    0,0,0,   -1*dst(1)*(src')]);
    
end

% Obtain SVD of A
[~,~,V] = svd(A);

% Find H using last column of V (the null space of A)
transform_mat = reshape(V(:,end), 3, 3)';
transform_mat = transform_mat./transform_mat(9);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Process Images
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Iterate through all of the images (only once if a single image!)
first_flag = 1;
for imnum = 1:numel(file_list)
      
    
    % Read in the target image
    if folder_flag
        input = imread([filepath file_list(imnum).name]);
    else
        input = imread(filepath);
    end
    
    % Guess at the bit depth of the image.
    if isa(input, 'uint8')
        tar_bitval = 8;
        input = im2double(input);
    elseif isa(input, 'uint16')
        tar_bitval = 16;
        input = im2double(input);
    else
        tar_bitval = ceil(log(double(max(input(:))))/log(2));
        input = input./(2^tar_bitval-1);
    end
    
    
    % Process the mosaic by splitting across color channels. Populate all other
    % pixels with -1 to avoid computational biasing.
    if mosaic_array == 0
        target = input;
    else
        tar_size = size(input);
        tar_height = tar_size(1);
        tar_width = tar_size(2);
        target = zeros([tar_size, 3])-1;
        target(1:2:tar_height, 1:2:tar_width, mosaic_array(1,1)) = input(1:2:tar_height, 1:2:tar_width);
        target(2:2:tar_height, 2:2:tar_width, mosaic_array(2,2)) = input(2:2:tar_height, 2:2:tar_width);
        target(1:2:tar_height, 2:2:tar_width, mosaic_array(1,2)) = input(1:2:tar_height, 2:2:tar_width);
        target(2:2:tar_height, 1:2:tar_width, mosaic_array(2,1)) = input(2:2:tar_height, 1:2:tar_width);
    end
    
    if draw_flag || first_flag
    
        close all;
        % Display the image, and read in 4 corners specified by user
        if peak_normalize == 0
            imshow(target);
        else
            imshow(target./max(target(:)));
        end
        hold on
        
    end
    
    
    % Iterate through each patch of the target
    for x_patch = 1:num_cols
        for y_patch = 1:num_rows
            
            % Define the current patch in the normalized coordinate system
            target_coords = [(x_patch-1 + crop_factor), (x_patch - crop_factor), (x_patch - crop_factor), (x_patch-1 + crop_factor); ...
                             (y_patch-1 + crop_factor), (y_patch-1 + crop_factor), (y_patch - crop_factor), (y_patch - crop_factor); ...
                                          1           ,                1         ,             1          ,             1         ];
            
            % Transform to target coordinate system
            img_coords = transform_mat*target_coords;
            
            % Convert from homogeneous coordinates to the images
            patch_coords = [img_coords(1,:)./img_coords(3,:); img_coords(2,:)./img_coords(3,:)];
            
            if draw_flag || first_flag
            
                % Draw the pretty red boxes
                plot([patch_coords(1,:), patch_coords(1,1)], [patch_coords(2,:), patch_coords(2,1)], '-r', 'LineWidth', 2);

                % Force a re-draw
                drawnow;
            
            end
            
            
            % Create a ROI mask
            mask = roipoly(target, patch_coords(1,:), patch_coords(2,:));
            
            for i = 1:num_channels
                
                % For each color channel, apply the mask
                patch = target(:,:,i);
                patch = patch(mask);
                patch = patch(patch ~= -1);
                
                % Mean the ROI and report
                output(y_patch, x_patch, i) = mean(patch(:));
                stdev(y_patch, x_patch, i) = std(patch(:));
                
            end
            
        end
    end
    
    if first_flag
        first_flag = 0;
        close all
    end
    
    if draw_flag
    
        figure;
    
        if peak_normalize == 0
            imshow(output, 'InitialMagnification', 'fit')
        else
            imshow(output./max(target(:)), 'InitialMagnification', 'fit');
        end

    else
        fprintf('File %d of %d.\n', imnum, numel(file_list));
    end
    
    %Create output values
    output = output.*(2^tar_bitval-1);
    stdev = stdev.*(2^tar_bitval-1);
    
    results(:,:,:,imnum) = output;
    standard_deviation(:,:,:,imnum) = stdev;
    
end

end

